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AN ADAPTIVE VARIABLE ORDER QUADRATURE STRATEGY 


PAUL HOUSTON AND THOMAS P. WIHLER 


Abstract. In this article we propose a new adaptive numerical quadrature procedure which 
includes both local subdivision of the integration domain, as well as local variation of the number 
of quadrature points employed on each subinterval. In this way we aim to account for local 
smoothness properties of the function to be integrated as effectively as possible, and thereby 
achieve highly accurate results in a very efficient manner. Indeed, this idea originates from 
so-called hp-version finite element methods which are known to deliver high-order convergence 
rates, even for nonsmooth functions. 


1. Introduction 


Numerical integration methods have witnessed a tremendous development over the last few 
decades; see, e.g., (2][3 151. In particular, adaptive quadrature rules have nowadays become an 
integral part of many scientific computing codes. Here, one of the first yet very successful ap¬ 
proaches is the application of adaptive Simpson integration or the more accurate Gauss-Kronrod 
procedures (see, e.g., (!!)• The key points in the design of these methods are, first of all, to keep 
the number of function evaluations low, and, secondly, to divide the domain of integration in such 
a way that the features of the integrand function are appropriately and effectively accounted for. 

The aim of the current article is to propose a complementary adaptive quadrature approach 
that is quite different from previous numerical integration schemes. In fact, our work is based on 
exploiting ideas from hp -type adaptive finite element methods (FEM); cf. [4)[§[T2)[T3)[20]. These 
schemes accommodate and combine both traditional low-order adaptive FEM and high-order (so- 
called spectral) methods within a single unified framework. Specifically, their goal is to generate 
discrete approximation spaces which allow for both adaptively refined subdomains, as well as 
locally varying approximation orders. In this way, the /ip-FEM methodology is able to resolve 
features of an underlying unknown analytical solution in a highly efficient manner. In fact, this 
approach has proved to be enormously successful in the context of numerically approximating so¬ 
lutions of differential equations, and has been shown to exhibit high-order algebraic or exponential 
convergence rates even in the presence of local singularities; cf. 9 17 18 


With this in mind, we adopt the /ip-adaptive finite element strategy for the purpose of intro¬ 
ducing a variable order adaptive quadrature framework. More precisely, we propose a procedure 
whereby the integration domain will be subdivided adaptively in combination with a local tuning 
of the number of quadrature points employed on each subinterval. To drive this refinement process, 
we employ a smoothness estimation technique from [6 22 (see also 12 for a related strategy), 


which was originally introduced in the context of /ip-adaptive FEM. Specifically, the smoothness 
test makes it possible to gain local information concerning the regularity of the integrand function, 
and thereby, to suitably subdivide the integration domain and select an appropriate number of 
quadrature points for each subinterval. By means of a series of numerical experiments we demon¬ 
strate that the proposed adaptive quadrature strategy is capable of generating highly accurate 
approximations at a very low computational cost. The main ideas on this new approach together 
with a view on practical aspects will be discussed in the subsequent section. 
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2. An ftp-T ype Quadrature Approach 


2.1. General Quadrature Rules. Typical quadrature rules for the approximation of an integral 


of a continuous function / : [—1,1] 


I ■= J f{x)dx 

l , take the form 

p 


i ~ Q p (f) — y ] Wp^kf(xpk)i 


( 2 . 1 ) 


( 2 . 2 ) 


fc=i 


where p > 1 is a (typically prescribed) integer number, and C [—1,1] and {w Pi fc}^ =1 C 

(0, 2] are appropriate quadrature points and weights, respectively. When dealing with a variable 
number p of quadrature points and weights, we can consider one-parameter families of quadrature 
rules (such as, for example, Gauss-type quadrature methods); here, for each p £ N, with p > 
Pmin, where p m \ n is a minimal number of points, there are (possibly non-hierarchical) families of 
quadrature points x p = {aip,fc}f = and weights w p = {w p ,fc}f =1 . 

On an arbitrary bounded interval [a, ft], a < ft, a corresponding integration formula can be 
obtained, for instance, by means of a simple affine scaling 


<l>[a,b] ■ [-1,1] -> [a, 6], 
with ft = ft — a > 0. Indeed, in this case 


^ x = <(>[a,b](£) = \hx+ *(a + ft), 


(2.3) 


fb h P 

/ f{x) ~ Q[a,b],p(f) • ^ ^ w p,k (/ ° ^[a.bl ) i.Xp,k\ 

1 k= 1 

where / : [a, b] —> R. is again continuous. As before, for any specific family of quadrature rules, 
the corresponding quadrature point families x p are obtained in a straightforward way by letting 
Xp = 4>[a,b](x p ) (with the understanding that 4>[ a ^] is extended componentwise to vectors). 

Furthermore, the above construction allows us to define composite quadrature rules, whereby 
the integral of / is approximated on a collection of n > 1 disjoint (open) subintervals {A'i}™ =1 
of [a, 6] with [a, b] = (J" = i K i: i.e., 




QKi,p(f\Ki)- 


In practical applications the subintervals are usually either of uniform size ( b ~ a )/n , for sufficiently 
large n, or alternatively, they are selected adaptively with the aim of resolving the relevant features 
of the given function /. 

2.2. The Basic Idea: ftp- Adaptivity. Adaptive quadrature rules usually generate a sequence 
of repeatedly bisected and possibly non-uniform subintervals {Ki}\ L 1 , n > 1, of the integration 
domain [a, ft] (i.e., each subinterval Ki may have a different length hi), with a prescribed and 
uniform number p of quadrature points on each subinterval. With the aim of providing highly 
accurate approximations with as little computational effort as possible, the novelty of the approach 
presented in this article is to design an adaptive quadrature procedure, which, in addition to 
subdividing the original interval [a, ft] into appropriate subintervals, is able to adjust the number 
of quadrature points pi individually within each subinterval K t in an effective way. We note that 
this idea originates from approximation theory HE! (see also and has been applied with huge 
success in the context of finite element methods for the numerical approximation of differential 
equations. Indeed, under certain conditions, the judicious combination of subinterval refinements 
(ft-refinement) and selection of local approximation orders (p- refinement), which results in the class 
of so-called ftp-finite element methods, is able to achieve high-order algebraic or exponential rates 
of convergence, even for solutions with local singularities; see, e.g. [18]. In an effort to automate 
the combined ft- and p-refinement process, a number of ftp-adaptive finite element approaches 
have been proposed in the literature; see, e.g, the survey article 14 and the references cited 


therein. In the current article, we pursue the smoothness estimation approach developed in 6.22 
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(cf. also 12 for a related methodology), and translate the idea into the context of adaptive variable 
order numerical quadrature. 

Starting from a subinterval Ki with pi quadrature points, we are given a current approxima¬ 
tion QKi,pi(f\Ki) of the subintegral 



f(x)dx 


QKi,Pi(f\Ki)- 


(2.4) 


Then, with the aim of improving the approximate value (/| k, ), in the sense of an /ip-adaptive 

finite element methodology in one-dimension, we propose two possible refinements of Kp. 

(i) /i-refinement: The subinterval Ki of length hi is bisected into two subintervals K} and Kf of 
equal size h */ 2 , and the number pi of quadrature points is either inherited to both subintervals 
or, in order to allow for derefinement with respect to the number of local quadrature points, 
reduced to pi — 1 points. In the latter case, we obtain a potentially improved approximation 


Q\,(f) = Q k \ ,max(l,pi — 1) (/) Q Kf ,max(l,pi —1) (/) 


(2.5) 


of fl2.4 >. 

(ii) p-refinement: The subinterval K, is retained, and the number pi of quadrature points pi is 
increased by 1, i.e., pi ■<— pi + 1. This yields an approximation 


Q p Ki (f) = Q Ki , Pi+ 1 (/). ( 2 . 6 ) 

In case that Pi = p max , where p max is a prescribed maximal number of quadrature points on 
each subinterval, we define 


Qx,(/) — QK hPi (f) + Qa' 2 , P ;(/)> (2-7) 

where Kj and Kf result from subdividing Ki as in ([i]). 

In order to determine which of the above refinements is more appropriate for a given subinterval Ki, 
we apply a smoothness estimation idea as outlined in the subsequent section. Once a decision 
between h- and p-refinement for Ki has been made, the procedure is repeated iteratively for any 
subintervals Ki for which Qxi,pi (f\K t ) and its refined value (resulting from the chosen refinement) 
differ by at least a prescribed tolerance tol > 0. 


2.3. Smoothness Estimation. The basic idea presented in the articles (6 12 22 is to estimate 
the regularity of a function to be approximated locally. Then, following along the lines of the 
/ip-approximation approach, if the function is found to be smooth, according to the underlying 
regularity estimation test, then a p-refinement is performed, otherwise an /i-refinement is employed. 
In |6|, the following smoothness indicator, for a (weakly) differentiable function / on an interval K :n 
has been introduced (cf. i Eq. (3)]): 


I L°°(Kj) 


?Kj [/] := hj 1/3 \\f\\ LHKj) + ||/'|| w 


V 2 


if f\K, ^ 0, 


if f\ Kj = 0. 


(F) 


The motivation behind this definition is the continuous Sobolev embedding W 1,2 (Kj) c —^ L°°(Kj), 
which implies that 

l|w|lioo(K- ) 

sup — zr . -—j-7- < 1; 

veHi( Kj ) h j /2 \\v\\ L2(Kj) + ^ hj 2 \\v'\\ L2{Kj) 

see 16, Proposition 1]. In particular, it follows that Tk, [/] < 1 in / is classified as being 
smooth on Kj if Tk 3 [/] > r, for a prescribed smoothness testing parameter 0 < r < 1, and 
nonsmooth otherwise. 

To begin, we first consider the special case when / is a polynomial of degree p ? > 1. Then, the 
derivative of order p 3 — 1 of / is a linear polynomial, and the evaluation of the smoothness 
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indicator [f < ' Pj 1 - ) ] from |f]) is simple to obtain. In fact, let us write /|^. in terms of a (finite) 
Legendre series, that is, 

Pj 

/k=Z>(W^})’ (2-8) 

1=0 

for coefficients ao, • ■ ■ ,a Pj £ R. Here, Li, l > 0, are the Legendre polynomials on [—1,1] (scaled 
such that Li{ 1) = 1 for all l > 0), and <f>Kj is the affine scaling of [—1,1] to K :) ; cf. ( |2.3| ). For / as 
in (2.8) it can be shown that 

1 + 


•Ac, 






1 + \Qj + ^2£; 


(2.9) 


■Pi 


where £ Pj = (2 pj — 1) | a Pj/a Pj -i| (provided that a Pj ._i ^ 0); see 16, Proposition 3]. In particular, 
this implies that 

1 ‘ ^ <Tk, |> 3 


V6 + 1 


-i) 


< 1 ; 


( 2 . 10 ) 


cf. (6} §2.2], 

In the context of the numerical integration rule ( |2.2| ) , the above methodology can be adopted as 
follows: suppose we are given pj > 2 quadrature points and weights, {Xp-kYk—i and {w p .k} p ^_ v 
respectively. Then, 

r h Pj 

/ f(x)dx « Q Kj , Pj (f\ Kj ) = -y o (j) Ko )(x Pjtk ). (2.11) 

JK i 1 k =1 

We denote the uniquely defined interpolating polynomial of / of degree pj — 1 at the given quad¬ 
rature points by 

Pj -1 

1=0 

Due to orthogonality of the Legendre polynomials, we note that 
21 + 1 


bi = 


hj 


IK, 


TLK jlPi -if(x)(Li o <t> K )(x) 6x, 1 = 0,... ,pj - 1 . 


We further assume that the quadrature rule under consideration is exact for all polynomials of 
degree up to 2 pj — 2. Thereby, 

Pi 


h = 


21 + 1 
2 

21 + 1 


55 w w>fc( n ^>w-i/) 0< M x Pj,k)Li(x pjt k) 


k= 1 
Pi 


5 . X Pi ,k ( f ° *pKj ) { x Pj (x P j t fc). 


k =1 


Consequently, we infer that 


£ K i,Pi~ 1 ' ~ (2 Pj 3) 


b Pi-i 


= (2 Pj - 1) 


Opi-2 

"Yhk =1 w Pj,k{f ° 4 > K j )(x pjt k)L pj -i(x pj ^k) 
1 w Pj,k(f ° < l ) Kj)(, x pj,k)^- l pj—2( x pj,k) 


( 2 . 12 ) 


and thus, in view of (2.9), we use the quantity 

1 + 


F Kj.pjif) : = 


_ _ £ (W 

1 + Ukj ,p :j -\ + ^ , Pj -i V + 1 


(2.13) 


cf. (2.10), to estimate the smoothness of f\x r Here, we emphasise that the computation of 1;k,, P j - 1 
does not require any additional function evaluations of / since the values (/ o (j> K,){ x p u k), k = 
1,... ,Pj, have already been determined in the application of the quadrature rule (2.11). 
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2.4. Adaptive Variable Order Procedure. Based on the above derivations, we now propose 
an hp -type adaptive quadrature method. To this end, we start by choosing a tolerance tol > 0, a 
smoothness parameter r £ y^/(V 6+i), l), and a maximal number p max > 2 of possible quadrature 
points on each subinterval. Furthermore, we define the interval K\ = [a, b], and a small number pi, 
2 < pi < p max , of quadrature points on K\. Moreover, we initialise the set of subintervals subs, the 
order vector p containing the number of quadrature points on each subinterval, and the unknown 
value Q of the integral as follows: 

subs = {K{\, p = {pi}, Q = 0. 


Then, the basic adaptive procedure is given as follows: 

1: while subs ^ 0 do 

2: [Ql, subs, p] = hprefine(/, subs, p,p max , r); 

3: Q = Q + Ql; 

4: end while 

5: Output Q. 

Here, hpref ine is a function, whose purpose is to identify those subintervals in subs, which need to 
be refined further for a sufficiently accurate approximation of the unknown integral. In addition, 
it outputs a set of subintervals (again denoted by subs), as well as an associated order vector 
(again denoted by p) which result from applying the most appropriate refinement, i.e., either h- 
or p-refinement as outlined in (i) and (ii) in Section 2.2 above, for each subinterval. Furthermore, 
hpref ine returns the sum Ql of all quadrature values corresponding to subintervals in the input 
set subs for which no further refinement is deemed necessary. The essential steps are summarised 
in Algorithm [l] 


Algorithm 1 Function [Q, subsnew, pnew] = hpref ine(/, subs, p,p max , r) 

1: Define subsnew = subs, and pnew = p. Set Q = 0. 

2: for each subinterval Kj £ subs do 
3: Evaluate the smoothness indicator F k from 

4: if (/) < r then 

5: Apply /i-refinement to Kj, i.e., bisect Kj into two subintervals of equal size and reduce 

the number of quadrature points to max(pj — 1,1) on both of them; 

6: Compute an improved approximation, denoted by Qx t , of Q k, , Pj (f \ x, ) using 

on Kj. 

7: else if F KjtPj (/) > r and Pj + 1 < p max then 

8: Apply p-refinement to Kj, i.e., increase the number of quadrature points to Pj + 1 

on Kj -, 

9: Compute an improved approximation, denoted by Qx, , of Q k, , p , if \ K, ) using 

on Kj. 

10: else if F Kj , Pj (/) > t and Pj + 1 > p max then 

11: Bisect Kj into two subintervals of equal size and retain the number of quadrature 

points pj on both of them; 

12: Compute an improved approximation, denoted by Qx t , of Q k, , p , if \ x, ) using 

on Kj. 

13: end if 

14: if \Qxj - Qx 0 , Pj {f IaTj)! is sufficiently small then 

15: Update Q = Q + Q k . ; 

16: Eliminate I\j from subsnew and the corresponding entry pj from pnew. 

17: else 

18: Replace Kj and pj in subsnew and pnew, respectively, by the corresponding h- or 

p-refined subintervals as determined above. 

19: end if 

20: end for 


(2.71 


(2.61 


(2.51 


(2.13). 
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2.5. Practical Aspects. In this section we discuss a number of practical issues involved in the 
implementation of the procedure described in Section [2.4| within a given computing environment. 


2.5.1. Gauss-Quadrature Rules. In principle, the adaptive procedure presented in Section 2.4 al¬ 
lows for any variable order family of quadrature rules. In our numerical experiments presented in 


Section 2.6 below, we propose the use of (families of) Gauss-type quadrature schemes. Although 
they might be criticised for their non-hierarchical structure, in the sense that they require more 
function evaluations in comparison to more traditional schemes (such as, for example, the adap¬ 
tive Simpson or fixed-order Gauss-Kronrod rules), our numerical results indicate that their high 
degree of accuracy may be exploited in a very efficient manner within the /ip-setting, particularly 
for smooth functions, with or without locally singular behaviour. Indeed, whilst non-hierarchical 
lower-order Gauss-type quadrature schemes might not be computationally competitive, it is a well- 
known feature of /ip-methods (see, e.g., [18| ) that their superiority becomes especially apparent 
on a variable, higher-order level. 

In the current article we employ Gauss-Legendre quadrature points and weights (with at 
least Pmin = 2 points and weights); these quantities can be precomputed up to any given or¬ 
der p max (in practice p max = 15 is usually more than sufficient) or even be generated on the spot 
in an efficient way (see, e.g., FTP1) if an upper bound p max cannot be fixed. In addition, we note 
that the Gauss-Legendre rule based on p points has a degree of exactness of 2p— 1, i.e., the smooth¬ 


ness indicators derived in Section 2.3 can be computed by means of the formula given in (2.12). 
For a given maximum number p max , we store the points and weights of the Gauss-Legendre rules 
(on the reference interval [—1,1]) with up to p max points in two p max x (p max — l)-matrices X 
and W, respectively; here, for parameters p = 2,... ,p max , the p-th columns of X and W are 
built from the points and weights of the corresponding p-point Gauss-Legendre quadrature rule, 
respectively (and complementing the remaining entries in all but the last column by zeros): 


(2.14) 


We note that, for other quadrature rules, the number of rows in the above matrices may be 
different. 



/^2,1 

X 3 ,l ' 

X v 1 \ 

maxi 1 ' 


(W 2 ,l 

w 3 ,1 • 

w v 1 \ 

y max ) 1 


^2,2 




W 2 ,2 



X = 


^3,3 


, W = 


U>3 ,3 



V 

0 

^Pmax iPmax ) 


\ 

0 

W V V / 

y max iy max / 


2.5.2. Vectorised Quadrature. Following the ideas of 19 we use a vectorised quadrature imple¬ 
mentation. This means that, instead of computing the integrals on the subintervals subs in 
Algorithm [l] one at a time, they are all computed at once. This can be accomplished by using fast 
vector- and matrix-operations, and by carrying out all necessary function evaluations in a single 
operation by computing the function to be integrated for a vector of input values. Specifically, we 
write the composite rule 


QK i ,p i {f\K i )= xlt, W PiAf ° 


6subs 


Ki(zsvibs 


k=l 


as a dot product of a weight vector w and a function vector f(x ); here, the former vector contains 
all (scaled) weights {kK w pi,k}i,k, and the latter vector represents the evaluation of the integrand 
function / on the vector x of all corresponding quadrature points {<j>Ki {x Pil k)}i.k appearing in 
the sum above. Evidently, these vectors can be built efficiently by extracting (and affinely map¬ 
ping and scaling) the corresponding rows from the matrices X and W in ( 2.14| ). We emphasise 
that applying vectorised quadrature crucially improves the performance of the overall adaptive 
procedure (provided that such a technology is available in a given computing environment). 
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2.5.3. Smoothness Estimators. As mentioned before, computing the smoothness indicators from 


(2.12) does not need any additional function evaluations of the integrand function /; they only 
require the values of the Legendre polynomials £ p _i and £ p _ 2 at the points {® P ,fe}fe =1) for p = 
2,..., Pmax• These quantities are again precomputable, and can be stored in two matrices 


and 


( L \( x 2 , i ) £2(23,1) 
£1(22,2) 

£2(23,3) 

0 

/£o( 22 ,l) £ 1 ( 23 , 1 ) 

£ 0 ( 22 , 2 ) 

£ 1 ( 23 , 3 ) 


^Pmax 1 (2p max ,1) ^ 


J-Jr) — 1 (l 
Kmax \ 


(2.15) 


Pmax iPmax 


)/ 


L n — 2 (*-£« 1 ) \ 

/-'max /-'max j ’ 


V 


0 


Ijir\ _ 0(0 

Kmax " \ 


(2.16) 


Pmax iPmax 


)/ 


Then, the sums in (2.12) are vectorised similarly as described above. In particular, the computation 


of the smoothness estimators can be undertaken with an almost negligible computational cost. 

2.5.4. Stopping Criterion. In order to implement the stopping-type criterion in line 14 of Algo¬ 
rithm [T] we exploit an idea that was proposed in the context of adaptive Simpson quadrature 
in 17 . More precisely, given a possibly rough approximation iguess ss J^ /( 2 ) da; of the unknown 
integral I from (2.1) (e.g., obtained from a Monte-Carlo calculation such that both the approx¬ 
imation and the exact value are of the same magnitude; cf. 7|), and a tolerance tol > 0, we 
redefine 

iguess = iguess * tol/eps; 

here, eps represents the smallest (positive) machine number in a given computing environment. 
Then, using the comparison operator ==, we accept the difference \Qkj ~ Qk.^p, (/|r, )| to be 
sufficiently small with respect to the given tolerance tol if the logical call 

iguess + | Q k - <5a,, p ,(/|a,)| == iguess; 


yields a true value. 


2.6. Numerical Examples. In order to test our approach, we consider a number of benchmark 
problems on the interval [0,1]. Specifically, the following functions will be studied: 


fx(x) = exp(x), 

/a(2) = V\ x ~ Vs|, 

fs{x) = sech(10(a; — V 5 )) 2 + sech(100(a; — 2 / 5 )) 4 

+ sech(1000(a: — 3 / 5 )) 6 + sech(1000(x — 4 /s)) 8 , 
fi(x) = cos(1000a;), 


h ( 2) 


0 if x < y3, 

1 if x > y 3 . 


Whilst the first function, / 1 , is analytic, the second function, / 2 , is smooth except at y.3 (see 
Figure [I] (top)). Furthermore, was proposed in [lO in the context of the chebfun package II 


this is a smooth function that exhibits several very thin spikes (see Figure [ 2 ] (top)). Moreover, 
is highly oscillating, and f 5 is an example of a discontinuous function. 
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# fct. calls 

hp-adapt. quad. 

# sing. fct. ev. 

cpu [sec] 

adapt. Simpson quad. 
# sing. fct. ev. 

h 

52 

9 

0.0031 

4,096 

h 

1,718 

65 

0.0224 

25,488 

h 

2,427 

33 

0.0144 

72,528 

h 

50,534 

35 

0.0180 

1,965,376 

fh 

1,273 

106 

0.0342 

784 


Table 1. Performance data for ftp-type adaptive quadrature. 


We perform our computations in MATLAI0 on a single 2.6GHz processor. The tolerance is set 
to tol = 0.3 x 10” 15 (which is close to machine precision in Matlab), the smoothness estimation 
parameter is prescribed as r = 0.6, and p m ax = 15. Within this setting, the adaptive proce¬ 
dure generates results that are accurate to machine precision, for all of the considered examples. 
In Table [T] for each of the functions /i,...,/ 5 above, we present the number of function calls 
(# fct. calls) in the vectorised quadrature implementation (counting a single application of the 
integrand function to a vector input as 1; cf. Section 2.5.2), as well as the number of single function 
evaluations (# sing. fct. ev.) taking into account the number of scalar entries of a vector input in 
each function call. The latter number is compared with the number of scalar function evaluations 
performed in a classical adaptive Simpson procedure as proposed in 7 (which is based on em¬ 
ploying the two end points as well as the midpoint on each subinterval, and reuses the former two 
points without recomputing). Except for the last function, /s, where a low-order quadrature rule 
is more effective, the remarkable efficiency of the proposed hp -type quadrature becomes clearly 
visible. This is confirmed with the expeditious cpu times (which do not include the computa¬ 
tion of the precomputable matrices X, W, L±, L 2 from (2.14), (2.15), and (2.16)) for each of the 
examples. 

In order to illustrate how the ftp-adaptive procedure performs, we depict the final ftp-mesh for fi 
and /*3 in Figure [l] (bottom) and Figure [2] (bottom), respectively. Here, along the horizontal axis 
we present the subintervals obtained as a result of the adaptive process, and on the vertical axis 
the number of quadrature points introduced on each subinterval is displayed. In both examples, we 
see that smooth regions in the underlying integrand are resolved by employing larger subintervals 
featuring a higher number of quadrature points, whereas close to singularities, the number of 
quadrature points is kept low on very small integration subdomains. It is noteworthy that this 
behaviour is well-known from ftp-finite element methods for differential equations, where high- 
order algebraic or even exponential convergence rates can be obtained by applying this type of 
ftp-refinement procedure; see 


18 for details. 


3. Conclusions 

In this article we proposed a new adaptive quadrature strategy, which features both local sub¬ 
division of the integration domain, as well as local variation of the number of quadrature points 
employed on each subinterval. Our approach is inspired by the ftp-adaptive finite element method¬ 
ology based on ftp-adaptive smoothness testing. In combination with a vectorised quadrature 
implementation, the proposed adaptive quadrature algorithm is able to deliver highly accurate 
results in a very efficient manner. Since our approach is closely related to the ftp-finite element 
technique, it can be extended to multiple dimensions, including, in particular, the application of 
anisotropic refinements of the underlying domain of integration, together with the exploitation of 
different numbers of quadrature points in each coordinate direction on each subinterval (based, 
for example, on anisotropic Sobolev embeddings as outlined in ® §3.1]). 


1 The MathWorks, Inc. 
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Figure 1 . Function / 2 : Graph (top) and hp -mesh (bottom). 
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5. R. DeVore and K. Scherer, Variable knot, variable degree spline approximation to x@, Quantitative approxima¬ 
tion (Proc. Internat. Sympos., Bonn, 1979), Academic Press, New York-London, 1980, pp. 121-131. 

6. T. Fankhauser, T. P. Wihler, and M. Wirz, The hp-adaptive FEM based on continuous Sobolev embeddings: 
isotropic refinements, Computers &; Mathematics with Applications. An International Journal 67 (2014), no. 4, 
854-868. 

7. W. Gander and W. Gautschi, Adaptive quadrature—revisited , BIT 40 (2000), no. 1, 84-101. 

8. A. Glaser, X. Liu, and V. Rokhlin, A fast algorithm for the calculation of the roots of special functions , SIAM 
Journal on Scientific Computing 29 (2007), no. 4, 1420-1438. 

9. W. Gui and I. Babuska, The h, p and h — p versions of the finite element method in one-dimension, parts 
I-III, Numer. Math. 49 (1986), no. 6, 577-683. 

10. N. Hale, Spike integral, 2010, http://www.chebfun.org/examples/quad/SpikeIntegral.html. 





ADAPTIVE VARIABLE ORDER QUADRATURE 


11 


11. N. Hale and L. N. Trefethen, Chebfun and numerical quadrature, Science China Mathematics 55 (2012), no. 9, 
1749-1760. 

12. P. Houston and E. Siili, A note on the design of hp-adaptive finite element methods for elliptic partial differ¬ 
ential equations, Comput. Methods Appl. Mech. Engrg. 194(2-5) (2005), 229-243. 

13. J. M. Melenk and B. I. Wohlmuth, On residual-based a posteriori error estimation in hp-FEM , Adv. Comp. 
Math. 15 (2001), 311-331. 

14. W. F. Mitchell and M. A. McClain, A comparison of hp-adaptive strategies for elliptic partial differential 
equations, ACM Transactions on Mathematical Software (TOMS) 41 (2014), 2:1—2:39. 

15. W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical recipes, third ed., Cambridge 
University Press, Cambridge, 2007, The art of scientific computing. 

16. K. Scherer, On optimal global error bounds obtained by scaled local error estimates, Numer. Math. 36 (1980/81), 
no. 2, 151-176. 

17. D. Schotzau, C. Schwab, and T. P. Wihler, hp-DGFEM for second-order mixed elliptic problems in polyhedra, 
Math. Comp, (in press). 

18. C. Schwab, p- and hp-FEM - Theory and application to solid and fluid mechanics, Oxford University Press, 
Oxford, 1998. 

19. L. F. Shampine, Vectorized adaptive quadrature in Matlab, J. Comput. Appl. Math. 211 (2008), no. 2, 131-140. 

20. P. Solin, K. Segeth, and I. Dolezel, Higher-order finite element methods, Studies in advanced mathematics, 
Chapman & Hall/CRC, Boca Raton, London, 2004. 

21. J. Waldvogel, Fast construction of the Fejer and Clenshaw-Curtis quadrature rules, BIT. Numerical Mathe¬ 
matics 46 (2006), no. 1, 195-202. 

22. T. P. Wihler, An hp-adaptive strategy based on continuous Sobolev embeddings, J. Comput. Appl. Math. 235 
(2011), 2731-2739. 

School of Mathematical Sciences, University of Nottingham, University Park, Nottingham, NG7 

2RD, UK 

E-mail address: Paul. Houston@nottingham .ac.uk 

Mathematisches Institut, Universitat Bern, Sidlerstrasse 5, CH-3012 Bern, Switzerland 

E-mail address : wihler@math. unibe. ch 



